On a two-dimensional model of generalized thermoelasticity with application

A 2D first order linear system of partial differential equations of plane strain thermoelasticity within the frame of extended thermodynamics is presented and analyzed. The system is composed of the equations of classical thermoelasticity in which displacements are replaced with velocities, complemented with Cattaneo evolution equation for heat flux. For a particular choice of the characteristic quantities and for positive thermal conductivity, it is shown that this system may be cast in a form that is symmetric t-hyperbolic without further recurrence to entropy principle. While hyperbolicity means a finite speed of propagation of heat waves, it is known that symmetric hyperbolic systems have the desirable property of well-posedness of Cauchy problems. A study of the characteristics of this system is carried out, and an energy integral is derived, that can be used to prove uniqueness of solution under some boundary conditions. A numerical application for a finite slab is considered and the numerical results are plotted and discussed. In particular, the wave propagation nature of the solution is put in evidence.

The hyperbolic systems of partial differential equations have been a subject of permanent interest, whether in the theory for their unique properties, or for their appicability to the phenomenon of wave propagation in several branches of Science and Technology. Problems of generalized thermoelasticity leading to heat wave propagation have provided great opportunity to investigate solutions of hyperbolic systems of partial differential equations under various initial and boundary conditions. Ruggeri 1 presented a survey on the relations between mathematical problems for quasi-linear hyperbolic systems and extended thermodynamics in continuum theories. It was shown that the system of balance laws can be symmetrized by use of an appropriate choice of the field variables. A purely thermal case was treated by Wilmański 2 , ch. 9] using an undetermined multipliers technique. Müller 3 investigated the symmetric hyperbolic systems of partial differential equations in extended thermodynamics. Selivanov and Selivanova 4 studied the computability properties of a system of symmetric hyperbolic equations by a difference scheme. Othman et al. 5 investigated generalized thermoelastic diffusion in a homogeneous, isotropic elastic half-space based on the Green-Naghdi theory. Abbas and Zenkour 6 used finite elements to investigate a two-dimensional problem under Green and Naghdi theory of thermoelasticity for a fiber-reinforcement anisotropic half-space subjected to a thermal boundary shock to assess the effects of initial stress and rotation. Cimmelli et al. 7 presented a review of the modern mathematical methods in the generalized thermodynamics of continuous media. He et al. 8 studied a two-dimensional generalized thermoelastic diffusion problem for a half-space. Mishra 9 solved a 2D-problem of heat transfer in a thin plate based on single-phase-lagging heat conduction by superposition technique and Fourier series expansion. Ghaleb et al. 10 presented a model of electrothermoelasticity within the theory of generalized thermodynamics. Abbas and Marin 11 studied a 2D-problem of generalized thermoelasticity for a half-space with surface laser heating. Rogolino et al. 12 derived two generalized heat-transport hyperbolic equations and studied thermodynamical compatibility for both. Jou 13 considered some fundamental aspects of non-equilibrium thermodynamics with a view on heat transport in nanosystems. Mahmoud et al. 14 investigated nonlinear heat wave propagation in rigid thermal conductors. Alzahrani et al. 15 studied a two-dimensional problem of a porous medium in extended thermodynamics using spectral method. Ahmed et al. [16][17][18] investigated two-dimensional problems of piezo-thermoelasticity within the dual-phase-lag model for layered media and a quarter-space. Solutions were presented based on normal modes technique, as well as numerical solutions by finite differences. Concerning the solution of systems of hyperbolic partial differential equations, Bonet et al. 19 proposed a novel computational framework for thermoelasticity based on the use of a 1. Equations of motion 2. Heat equation 3. Cattaneo-Vernotte relations. These evolution laws for the heat flux components replace the classical Fourier law for heat conduction. Divided throughout by thermal relaxation times τ 1 and τ 2 respectively, these relations read: 4. The generalized Hooke's law differentiated w.r.t. time, thus allowing to eliminate the mechanical displacement components in favour of the velocities: The problem thus reduces to the solution of eight basic partial differential equations of the first order. These equations involve eight unknowns: two velocity components, three identically non vanishing stress components, temperature and two heat flux components. Having resolved this system of equations, one can then find the mechanical displacement components by quadrature from the relations: It is not difficult to anticipate a time damping of the heat wave, due to the terms involving n 11 and n 22 in the Cattaneo evolution equations. In order to relate the coefficients m 11 , m 22 , n 11 and n 22 in Eqs. (4 and 5) to known physical quantities, it is sufficient to compare these two equations to the well-known standard forms of the evolution equations for the heat flux (C.f. 10,21 ): where τ 1 , τ 2 are the thermal relaxation times for the x-and y-directions respectively, and k 11 , k 22 are the thermal conductivities in these two directions. Direct comparison of these last two equations with (4) and (5) reveals that The above governing equations need to be cast in a convenient form for later work. For this, Eqs. (6 and 7) will be replaced by their symmetrized forms: where It may be easily verified that Excluding Eqs. (9 and 10) which may be considered independently in a later stage as explained above, the system of eight basic equations is written in matrix form as: Matrices A, B and C have dimension (8 × 8) and are given as: In the case of multiple roots, the system still preserves all the main properties of hyperbolic systems.

Definition 3
The system of Eq. (15) is said to be symmetric, t-hyperbolic if the matrices A, B, C are symmetric and, moreover, matrix A is positive definite.
It is well-known that for symmetric, t-hyperbolic systems one can deduce the so-called energy integral, which represents a powerful tool to prove a theorem on the uniqueness of solution.

Reduction to symmetric t-hyperbolic form
In what follows, we investigate the possibility of diagonalizing the matrix τ A + ξ B + ηC for all values of (ξ , η) and τ > 0 , in order to simplify the considered system of equations and reveal its nature. As a first step, the basic equations will now be reformulated in dimensionless form. To this end, introduce the following set of dimensionless variables for length, time, temperature heat flux, velocity and stress: The characteristic quantities θ 0 , L 0 , T 0 , Q 0 , c are kept undetermined at this stage, but will be defined later on step by step so as to achieve some properties and simplifications. The two characteristic quantities L 0 and T 0 are related to each other by: It may be noticed that the characteristic quantities θ 0 and Q 0 are independent of each other according to the basic assumption of extended thermodynamics adopted in the present model (see 10 , for example). This fact provides some flexibility in attaining the required final form of the equations.
After removing the tilde, the system of eight basic equations in dimensionless form reads:  The matrices of the system now read: Matrices B and C can be made symmetric by suitable choices of some parameters, taking in consideration that some characteristic quantities of the system are still kept arbitrary. First, let us choose the characteristic velocity to be equal to the speed of propagation of the purely elastic transverse wave in the linear approximation: Next, assume the material parameters are such that the following relation holds: which amounts to defining the reference temperature as: The characteristic time T 0 will be determined depending on the thermal relaxation time in the numerical example treated below. The final form of the basic system of eight equations is: where for the transversely isotropic case under consideration. By elimination, one easily deduces from the above equations the following expression for the dimensionless speed of the heat wave:

Characteristics
In order to find the characteristics of the system of partial differential equations under consideration, note that the symmetric matrices A, B and C are now given by: Introduce the matrix T whose columns are the eigenvectors of matrix A, arranged as shown above: The similarity transformation with matrix T: X ′ = T T XT applied to matrices A, B and C yields the following:  Then the two symmetric matrices A ′ and ξ B ′ + ηC ′ are simultaneously diagonalizable by means of the transformation S = A ′−1/2 O for all real values of (ξ , η) . Moreover, the diagonal form of the first matrix is the unit matrix: The diagonal matrix D may be obtained from the eigenvalues of matrix W. The characteristic equation of this matrix is: The expressions for the coefficients are as follows: from which it is seen that they are all positive. For sufficiently small absolute values of ξ and η , the value of α 1 is much larger than the values of the other two coefficients α 2 and α 3 . Consider the related cubic equation: One easily verifies that The local extrema of function f are located at: The discriminant for this cubic equation is: It is easy to see that sufficient conditions to get a curve shape for the cubic polynomial as in fig.1, i.e. for the eigenvalues to be real, are: For the used values of the material parameters, these two conditions are likely to be satisfied for all values of ξ and η . The fulfillment of the two conditions was achieved numerically as illustrated in Figs. 2 and 3. At the particular point (ξ = 0, η = 0) , all eigenvalues are coincident and equal to zero. Thus the Eq. (35) has eight real roots: A double root equal to zero, three positive roots, and three negative roots, the roots being symmetrically positioned with respect to the origin. One is finally led to the following:

The energy integral
It is well-known (c.f. 20 ) that an energy integral can be derived for symmetric t-hyperbolic systems of partial differential equations. If V denotes a region of the t, x, y -space included in the domain of definition of the solution, bounded by a surface S, then the energy integral for the transformed system reads 20 , p. 123]:   where n = (ξ , η) and we have introduced σ ′′ with components: It is known that the characteristic surfaces for the wave equation are parts of conuses with axis parallel to the t-axis, so that τ > 0 . The boundary conditions are taken in the form: yielding It is well-known that this last inequality leads to uniqueness of solution. The first two boundary conditions in (43) mean impermeability and thermal insulation of the boundary. The third one englobes many cases, among which complete fixing of the boundaries.

Numerical application
As a confirmation of the existence of solutions to the investigated system of equations, we have considered a Cauchy problem for rectangular slab under specified boundary conditions. The solution includes the propagation of three types of waves, the heat wave, the transversal and the longitudinal thermomechanical coupled waves. Setting τ = τ 1 = τ 2 the thermal relaxation time and k = k 11 = k 22 the coefficient of heat conduction, one has to make assumptions concerning the orders of magnitude of constants appearing in (4) and (5): Thus the dimensionless parameters M and n appearing in the dimensionless Cattaneo Eqs. (32 and 33) satisfy the rules For the example treated below, one has M ∼ 1 and n ∼ 1 , hence With this in mind, let us consider the case of a transversely isotropic material having tentative values of the different geometrical and material parameters shown in Table 1. Consultation of the available tables of material (43) v.n = q.n = σ ′′ v .n = 0,  www.nature.com/scientificreports/ coefficients for metals and alloys reveals that the used values of mass density and specific heat capacity as displayed in Table 1 are within the normal range, while the thermal conductivity lies in the higher range. Corresponding to these values: Motion is induced by a thermal boundary regime. For this, consider a rectangular domain occupied by a thermoelastic material, with dimensions a and b and sides labeled S 1 , S 2 , S 3 and S 4 as shown in Fig. 4. The origin of the plane Cartesian coordinate system is chosen at the left lower corner of the rectangle, with x-axis along the side S 4 . Initially, the medium is at rest at zero temperature (measured from an initial ambient temperature θ 0 ) and heat flux, and in a stress-free state. The chosen boundary conditions are of mixed type. Mechanically, sides S 1 , S 2 and S 4 are traction-free, side S 3 is completely fixed so that to suppress rigid body motion of the slab. As concerns the thermal conditions, the normal heat flux is taken to vanish on boundaries S 2 , S 3 and S 4 , together with Robin thermal condition. Side S 1 has a prescribed heating regime that generates the motion. It should be remembered that temperature and heat flux are independent thermodynamical quantities. For definiteness: (1) At the side S 1 (x = 0, 0 ≤ y ≤ b): and The graphical representation of this function for b = 0.4 is illustrated in Fig. 5.
(2) At the side S 2 (y = b, 0 ≤ x ≤ a): together with Robin thermal condition and Robin thermal condition and Robin thermal condition Here, Bi is a dimensionless Biot number, with value taken as Bi = 0.0144 . All the remaining functions not specified in the above boundary conditions are set to zero on the boundaries. Following similar guidelines as in 22,23 , COMSOL Multiphysics is used to solve the considered boundary-value problem for Eqs. (26-33) under the boundary conditions (46-52) and zero initial conditions. For definiteness, we mainly consider the case with M = 3.0, n = 1.0 . The value of M, however, was varied to take on other values in some plots. The produced particular solution has been obtained by the method of finite elements embedded in the above-mentioned software. The mesh could be refined and adjusted for best results.
σ S 3 xx a, y, t = σ S 3 xy a, y, t = q S 3 x a, y, t = u x a, y, t = u y a, y, t = 0 www.nature.com/scientificreports/ From symmetry of the considered problem about the median line y = b/2 , it follows that the quantities u y , v y , q y , σ xy must vanish everywhere in the slab. The numerical computations have confirmed this fact . It thus follows that the displacements, velocities and heat flux take place along the direction perpendicular to the heated face of the slab (the x-direction), and that shear stress vanishes inside the slab. Figures 6,7,8,9 ,10 and 11 represent topviews of the 3-D distributions at three consecutive time moments of the main physical quantities induced inside the material due to applied boundary conditions. The calculations are performed for a dimensionless heat wave speed M = 3 . This wave thus travels three times faster than the transversal coupled thermoelastic wave whose dimensionless speed is c = 1.0 . For the chosen values of the material constants, the longitudinal coupled thermoelastic wave travels at speed ≃ 2.944 . In all these figures, wave propagation phenomenon is clearly illustrated for both the mechanical and the thermal variables, with oscillations between positive and negative values near the heated side of the slab, and amplitudes monotonically decreasing to zero with distance. For the considered time values, the wave fronts are clearly noticed and the propagating disturbances have not yet reached the far end of the slab and therefore no reflected waves are expected. Larger time values could not be considered due to lack of stability of the computational scheme.
The first remark concerns the distribution of temperature in the slab. Although the motion was initiated by heating of the left boundary for a finite lapse of time, it is noticed in Fig. 6 that temperature assumes negative values near this boundary for the considered values of time. This seemingly paradoxical situation may be explained by the fact that part of the thermal energy supplied to the medium at the boundary is spent to generate mechanical wave propagation and, therefore, cooling may take place in some parts of the medium. Again, as shown by    24 for the case of a rigid thermal conductor, negative temperatures can result from unphysical values assigned to the thermal relaxation time.
As noted above, heat flux and temperature are independent thermodynamical variables. Heat flux is initially generated by the thermal boundary condition. In subsequent time moments, it propagates as a wave, with speed M, independently of temperature as shown in Fig. 7.
Still at the heated end, the displacement and the velocity components u x and v x in Figs. 8 and 9 take on negative values in the initial phase of the motion, as the medium expands under heating. Both quantities are seen to tend to zero as the corner points are approached.

Conclusions
We have investigated a two-dimensional system of first-order, partial differential equations of classical thermoelasticity, supplemented with Cattaneo-type evolution equation for the heat flux to replace Fourier law for heat conduction. The model includes only one thermal relaxation time for each component of the heat flux vector. It differs from other well-known models of extended thermodynamics, e.g. Lord and Shulman, Green and Naghdi, two-temperature model, dual-phase-lag model, which all start with different assumptions and yield heat wave propagation (C.f. 25 ). Under the assumption that thermal conductivity must be positive, it turns out that the present system is reducible to symmetric t-hyperbolic form by special choices of some characteristic quantities. This result is to be contrasted with that expressed by Müller 3 , stating that the set of quasi-linear first order balance equations of extended thermodynamics may be written in symmetric hyperbolic form by a suitable choice of fields, provided constitutive functions are subject to the requirements of the entropy principle. The characteristics of the system were studied and an energy integral was obtained which leads to uniqueness of solution under proper boundary conditions. Thus the present system of equations behaves well in the sense of well-posedness of Cauchy problems, and is therefore valid for the description of heat wave propagation. To confirm the existence of solutions, a numerical experiment for a finite slab was carried out using a finite element scheme built in COMSOL Multiphysics and with tentative values of the different material parameters to produce a particular solution to the problem under specified boundary conditions and zero initial conditions. The numerical results could be obtained only for sufficiently small time values due to difficulties related to the stability of the numerical scheme. The presented three-dimensional plots clearly show the phenomenon of thermal and thermomechanical wave propagation. It was noticed that temperature attains negative values at certain locations inside the slab, although the left boundary is subjected to heating for a certain lapse of time. This is due to the fact that part of the supplied thermal energy is spent on generating the motion. Such phenomenon should not appear in rigid thermal conductors once the value of the thermal relaxation time has been properly assigned (C.f. 24 ). The presented figures show the wave front and the damping of the solution as disturbances progress along the slab. Verification of the satisfaction of the boundary conditions was carried out.

Data availability
All data generated or analyzed during this study are included in this published article